import numpy as np
import matplotlib.pyplot as plt

plt.ion()

def sigt(j, E):
    return np.sin(j * np.log(E)) + 1.1

def sigs(j, E):
    return sigt(j, E) * 0.4 * (np.sin(4*j*np.log(E)) + 1)

E = np.logspace(-6, 0, 1000)

plt.semilogx(E, sigt(1, E), label='sigma t')
plt.semilogx(E, sigs(1, E), label='sigma s')
plt.semilogx(E, sigt(1, E) - sigs(1, E), label='sigma a')
plt.semilogx(E, sigs(1, E) / sigt(1, E), label='sigma s / sigma t')

plt.ylim(0,3)
plt.legend()
plt.savefig('xs-plot.png')

Sigt, Sigs = 0., 0.
for j in range(1,101):
    atom_density = 0.02 * np.exp(-0.001 * j)
    Sigt += sigt(j, E) * atom_density
    Sigs += sigs(j, E) * atom_density
plt.figure()
plt.semilogx(E, Sigt, label='sigma t')
plt.semilogx(E, Sigs, label='sigma s')
plt.semilogx(E, Sigt - Sigs, label='sigma a')
plt.semilogx(E, Sigs / Sigt, label='sigma s / sigma t')
plt.legend()
plt.savefig('macro-plot.png')
